/*
Copyright (c) [2019年5月1日] [吴超]
[MBT_studio] is licensed under Mulan PSL v2.
You can use this software according to the terms and conditions of the Mulan PSL v2.
You may obtain a copy of Mulan PSL v2 at:
		 http://license.coscl.org.cn/MulanPSL2
THIS SOFTWARE IS PROVIDED ON AN "AS IS" BASIS, WITHOUT WARRANTIES OF ANY KIND,
EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO NON-INFRINGEMENT,
MERCHANTABILITY OR FIT FOR A PARTICULAR PURPOSE.
See the Mulan PSL v2 for more details.
*/
#include "stdafx.h"
#include "GPU解算/PBF_Simulation.h"
#include "PBD解算器.h"


#include <ppl.h>
#include <几何图形/图形相交.h>
//#include <device_launch_parameters.h>
//#include <atomic_functions.h>
#include "graph_Adjacency.h"
#include "GPU解算/PBF_Simulation.h"
#include "GPU解算/PBF刚体解算.h"
#include "core/体素坐标生成.h"


#include "底层绘图/底层绘图框架.h"
#include <线性代数/几何计算.h>

#include <SIMD/向量运算.h>
#include <SIMD/CPU_向量计算.h>


using namespace std;



S_PBF解算器::S_PBF解算器() : m_邻接(E_并行设备::e_CUDA并行) {
	m_IdType = E_物理引擎::E_物理引擎_PBF;
	m_是否重构 = true;
	m_使用GPU计算 = true;

	
	m_物体添加有效粒子数量 = 0;
	m_最大迭代次数 = 3;
	m_更新步长 = 0.001;

	//m_邻接.m_Data.m_Num = 0;
	//m_邻接.m_Data.m_域属性.m_单元属性.m_分辨率 = { 1,1,1};
	//auto ctx	= f_创建协处理器设备(0, E_设备类型::e_计算设备);
	//m_计算设备	= f_创建计算环境(ctx);

	//m_外力计算 = new S_PBF外力计算(ctx);
	//m_解算数据 = new S_坐标数据(ctx);
	//m_邻接计算 = new S_邻接索引计算(ctx);
	m_飞溅粒子数量 = 0;
}

void S_PBF解算器::f_Update() {
	/*if (m_是否重构) {
		f_重构物理体();
		m_是否重构 = false;
	}*/
	/*if (m_是否重构) {
		f_重新计算各个物体的粒子索引偏移();

		m_粒子数据.m_流体粒子总数 = (m_物体添加有效粒子数量 / m_粒子数据.m_GPU数据.m_线程块大小) * m_粒子数据.m_GPU数据.m_线程块大小 + m_粒子数据.m_GPU数据.m_线程块大小;
		m_粒子数据.f_resize(m_粒子数据.m_流体粒子总数, true);


		vec3 空坐标 = m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移 + m_邻接.m_Data.m_域属性.m_单元属性.m_单元大小 * 2;
		vec3* p = m_粒子数据.m_位置.data();
		float32* w = m_粒子数据.m_质量.data();
		uint32* pID = m_粒子数据.m_物理体ID.data();
		uint8* pType = m_粒子数据.m_物理体类型.data();
		

		uint64 流体粒子数粒子数 = 0;
		uint64 offset = 0;
		uint32 物体ID = 0;
		uint64 num = 0;
		for (auto& body : m_流体) {
			S_PBF流体* b = (S_PBF流体*)body->m_Data;
			S_体素* v = b->m_体素;

			vec3* 体素坐标 = v->m_坐标.data();
			num = b->m_Count - b->m_Offset;
			流体粒子数粒子数 += num;

#pragma omp parallel for
			for (int64 i = 0; i < num; ++i) {
				uint64 id = b->m_Offset + i;

				p[id] = 体素坐标[i];
				w[id] = 0.8;

				pID[id] = 物体ID;
				pType[id] = body->m_Type;
			}

			++物体ID;
			//offset += v->m_坐标.size();
		}

		for (auto& body : m_流入) {
			S_PBF粒子发射* b = (S_PBF粒子发射*)body->m_Data;
			S_体素* v = b->m_体素;

			num = b->m_Count - b->m_Offset;
			流体粒子数粒子数 += num;

#pragma omp parallel for
			for (int64 i = 0; i < num; ++i) {
				uint64 id = b->m_Offset + i;

				p[id] = 空坐标;
				w[id] = 0.8;

				pID[id] = 物体ID;
				pType[id] = E_物理体类型::E_物理体_流入;
			}

			++物体ID;
		}


		uint64 刚体粒子数 = 0;
		std::vector<uint32> 刚体区间起始;
		std::vector<uint16> 刚体粒子数量;
		for (auto& body : m_刚体) {
			S_PBF刚体* b = (S_PBF刚体*)body->m_Data;
			S_体素* v = b->m_体素;

			vec3* 体素坐标 = v->m_坐标.data();
			num = b->m_Count - b->m_Offset;
			刚体粒子数 += num;

#pragma omp parallel for
			for (int64 i = 0; i < num; ++i) {
				uint64 id = b->m_Offset + i;
				p[id] = 体素坐标[i];
				w[id] = 0.8;

				pID[id] = 物体ID;
				pType[id] = body->m_Type;
			}

			++物体ID;

			//f_初始化粒子刚体变换(b->m_旋转矩阵, b->m_中心坐标, &(p[b->m_区间.x]), &(w[b->m_区间.x]), num);
			刚体区间起始.push_back(b->m_Offset);
			刚体粒子数量.push_back(num);
			//offset += v->m_坐标.size();
		}



		//处理多出部分
#pragma omp parallel for
		for (int64 i = m_物体添加有效粒子数量; i < m_粒子数据.m_流体粒子总数; ++i) {
			w[i] = 1.0;
			p[i] = m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移 + m_邻接.m_Data.m_域属性.m_单元属性.m_单元大小*2;
			//p[i] = {};
			pID[i] = 物体ID;
			pType[i] = E_物理体类型::E_物理体_静态;
		}


		
		
		if (m_使用GPU计算) {
			m_粒子数据.f_初始数据到GPU();
			f_CUDA_update刚体区间(m_粒子数据.m_GPU数据, 刚体区间起始.size(), 刚体区间起始.data(), 刚体粒子数量.data(), 0);
			//f_CUDA_拷贝刚体粒子初始坐标(m_粒子数据.m_GPU数据, 刚体粒子数, &(p[流体粒子数粒子数]));
		}


		m_是否重构 = false;
	}
	
	//构建静态网格碰撞划分
	if (m_使用GPU计算) {
		//std::cout << "f_CUDA_计算施加外力和速度叠加" << std::endl;
		//std::cout << "f_构建粒子GPU邻接单元" << std::endl;
		for (auto& body : m_流入) {
			f_GPU_添加外力(m_粒子数据.m_GPU数据, *((S_PBF粒子发射*)body->m_Data));
		}


		f_CUDA_计算施加外力和速度叠加(
			m_粒子数据.m_GPU数据,
			m_粒子数据.m_重力,
			{ 0.0001, 0.005 },
			m_更新步长,
			m_粒子数据.m_流体粒子总数,
			//邻接索引,
			m_邻接.m_Data.m_域属性.m_单元属性,
			m_使用GPU计算);

		
		for (auto& body : m_流入) {
			f_GPU_置粒子位置(m_粒子数据.m_GPU数据, *((S_PBF粒子发射*)body->m_Data));
		}


		f_GPU_坐标钳到域内(m_粒子数据.m_GPU数据, m_邻接.m_Data.m_域属性.m_单元属性, m_粒子数据.m_流体粒子总数);

		m_邻接.f_构建粒子GPU邻接单元(m_粒子数据.m_GPU数据.m_cu位置, m_粒子数据.m_流体粒子总数, m_粒子数据.m_GPU数据.m_线程块大小);
		//m_邻接.f_toCpu();

		


		

		S_cuAdjacencyIndex 邻接索引 = f_getGPU数据();
		for (uint32 i = 0; i < 3; ++i) {
			f_CUDA_计算投影约束(m_粒子数据.m_GPU数据, 邻接索引, m_邻接.m_Data.m_域属性.m_单元属性, m_粒子数据.m_流体粒子总数);
		}
		
		//处理刚体柔体
		f_CUDA_update刚体粒子(m_粒子数据.m_GPU数据, m_刚体.size());


		m_粒子数据.f_数据到CPU();

		f_CUDA_更新粒子速度(m_粒子数据.m_GPU数据, 邻接索引, m_邻接.m_Data.m_域属性.m_单元属性, m_粒子数据.m_流体粒子总数);

		
		f_计算静态网格碰撞();
		m_粒子数据.f_数据到GPU();

		//std::cout << "m_邻接.f_toCpu" << std::endl;
	}
	else {
		f_施加外力();

		m_邻接.f_构建粒子CUP邻接单元(m_粒子数据.m_位置.data(), m_粒子数据.m_流体粒子总数, m_粒子数据.m_GPU数据.m_线程块大小);

		f_投影约束();
		
		//更新速度
		f_更新速度矢量();

		//计算粘度
		f_解算粘度();

		f_计算静态网格碰撞();
	}*/
	
	return;
}

void S_PBF解算器::f_应用物理变换() {
	
}

//void S_PDB解算器::f_取更新物体(std::vector<std::vector<Mat44f>>& mat, std::vector<std::vector<void*>>& userData) {
//}

void S_PBF解算器::f_取碰撞物体(std::vector<S_碰撞信息>& objs, uint32 scenekey) {

}

void S_PBF解算器::f_添加物理体(S_物理体* body, uint32 scenekey) {
	uint64 begin = m_物体添加有效粒子数量;

	switch (body->m_Type) {
	case E_物理体类型::E_物理体_流体: {
		S_PBF流体* b = (S_PBF流体*)body->m_Data;

		b->m_Offset = begin;
		begin += b->m_体素->m_坐标.size();
		b->m_Count = begin;

	} break;

	case E_物理体类型::E_物理体_刚体: {
		S_PBF刚体* b = (S_PBF刚体*)body->m_Data;

		b->m_Offset = begin;
		begin += b->m_体素->m_坐标.size();
		b->m_Count = begin;
		break;
	}

	case E_物理体类型::E_物理体_流入: {
		S_PBF粒子发射* b = (S_PBF粒子发射*)body->m_Data;

		b->m_Offset = begin;
		begin += b->m_数量;
		b->m_Count = begin;
		break;
	}

	case E_物理体类型::E_物理体_静态:
		//m_静态Body.push_back(dynamic_cast<S_PBF物理体*>(body));
		break;
	default:
		break;
	}

	m_物体添加有效粒子数量 = begin;
	//begin = (begin / m_粒子数据.m_GPU数据.m_线程块大小) * m_粒子数据.m_GPU数据.m_线程块大小 + m_粒子数据.m_GPU数据.m_线程块大小;
	//m_粒子数据.f_resize(begin, true);

}

void S_PBF解算器::f_移除物理体(S_物理体* body, uint32 scenekey) {
	switch (body->m_Type) {
	case E_物理体类型::E_物理体_刚体:
		m_刚体.erase(body);

		delete ((S_PBF刚体*)(body->m_Data))->m_体素;
		delete (S_PBF刚体*)(body->m_Data);
		break;
	case E_物理体类型::E_物理体_流体:
		m_流体.erase(body);

		delete ((S_PBF流体*)(body->m_Data))->m_体素;
		delete (S_PBF流体*)(body->m_Data);
		break;

	case E_物理体类型::E_物理体_流入:
		m_流入.erase(body);

		delete ((S_PBF流体*)(body->m_Data))->m_体素;
		delete (S_PBF流体*)(body->m_Data);
		break;
	}

	m_是否重构 = true;
}

void S_PBF解算器::f_添加流体(_S_PBF流体* body) {

}

void S_PBF解算器::f_添加物理体碰撞回调(S_物理体* body, uint32 scenekey) {

}






void S_PBF解算器::f_销毁物理体(S_物理体* body) {
	switch (body->m_Type) {
	case E_物理体类型::E_物理体_刚体:
		m_刚体.erase(body);

		delete ((S_PBF刚体*)(body->m_Data))->m_体素;
		delete (S_PBF刚体*)(body->m_Data);
		break;
	case E_物理体类型::E_物理体_流体:
		m_流体.erase(body);

		delete ((S_PBF流体*)(body->m_Data))->m_体素;
		delete (S_PBF流体*)(body->m_Data);
		break;

	case E_物理体类型::E_物理体_流入:
		m_流入.erase(body);

		delete ((S_PBF流体*)(body->m_Data))->m_体素;
		delete (S_PBF流体*)(body->m_Data);
		break;
	}
	
}

void S_PBF解算器::f_重新计算各个物体的粒子索引偏移() {
	uint64 begin = 0;

	for (auto& body : m_流体) {
		S_PBF流体* b = (S_PBF流体*)body->m_Data;
		
		b->m_Offset = begin;
		begin += b->m_体素->m_坐标.size();
		b->m_Count = begin;
	}

	for (auto& body : m_流入) {
		S_PBF粒子发射* b = (S_PBF粒子发射*)body->m_Data;

		b->m_Offset = begin;
		begin += b->m_数量;
		b->m_Count = begin;
	}

	for (auto& body : m_刚体) {
		S_PBF刚体* b = (S_PBF刚体*)body->m_Data;

		b->m_Offset = begin;
		begin += b->m_体素->m_坐标.size();
		b->m_Count = begin;
	}

	m_物体添加有效粒子数量 = begin + m_粒子数据.m_飞溅粒子数量;
}













void S_PBF解算器::f_施加外力() {
	
	vec3* 加速度 = m_粒子数据.m_加速度.data();
	vec3* 外力 = m_粒子数据.m_外力.data();
	float32* 重量 = m_粒子数据.m_质量.data();
	auto* delta = m_粒子数据.m_位置量变化.data();


#pragma omp parallel for
	for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
		delta[i] = { 0,0,0 };

		//流体部分
		if (重量[i] != 0.0f) {
			外力[i] = m_粒子数据.m_重力;
			
			//其它力施加
			加速度[i] = m_粒子数据.m_重力;
		}
		//刚体部分
	}


	vec3* 位置 = m_粒子数据.m_位置.data();
	vec3* 速度 = m_粒子数据.m_速度.data();
	vec3* old位置 = m_粒子数据.m_上一帧位置.data();
	


	//计算更新步长
	f_更新解算步长(0.0001, 0.005);



//#pragma omp parallel for
	/*for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
		old位置[i] = 位置[i];

		if (重量[i] != 0.0f) {
			//速度[i] = vec_add(速度[i], 加速度[i] * m_更新步长);
			//位置[i] = vec_add(位置[i], 速度[i] * m_更新步长);
			速度[i] = vec_add(速度[i], 加速度[i] * 0.015);
			位置[i] = vec_add(位置[i], 速度[i] * 0.015);

			位置[i] = vec_Max(位置[i], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移);
			位置[i] = vec_Min(位置[i], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移 + m_邻接.m_Data.m_域属性.m_单元属性.m_世界大小);
		}
	}*/


}



void S_PBF解算器::f_投影约束() {
	vec3* 位置 = m_粒子数据.m_位置.data();
	float32* 密度 = m_粒子数据.m_密度.data();
	float32* 质量 = m_粒子数据.m_质量.data();
	float32* 拉格朗日 = m_粒子数据.m_Lambda.data();
	vec3* 位置量变化 = m_粒子数据.m_位置量变化.data();


	for (uint8 迭代 = 0; 迭代 < m_最大迭代次数; ++迭代) {
		cout << "" << endl;
		this;
#pragma omp parallel for
		for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
			if (i == 1891) {
				float32 a = 0;
				a += 1;
				位置[i];
			}
			f_解算密度(i, 位置, 质量, m_粒子数据.m_密度0, 密度);
			拉格朗日[i] = f_计算拉格朗日乘数(i, 位置, 质量, 密度[i], m_粒子数据.m_密度0);

			if (i == 1891) {
				cout << "拉格朗日[i] :" << 拉格朗日[i] << endl;
				cout << "拉格朗日[i] :" << 密度[i] << endl;
			}
		}


		//解算密度约束
#pragma omp parallel for
		for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
			vec3 c = f_解算密度约束(i, 位置, 质量, 拉格朗日, m_粒子数据.m_密度0);
			if (i == 1891) {
				cout << "C x:" << c.x << " y:" << c.y << " z:" << c.z << endl;
			}
			//vec_add(&位置量变化[i], c);
			位置量变化[i] = c;
		}
		

#pragma omp parallel for
		for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
			
			if (i == 1891) {
				cout << "位置量变化[i] x:" << 位置量变化[i].x << " y:" << 位置量变化[i].y << " z:" << 位置量变化[i].z << endl;
			}
			vec_add(&位置[i], 位置量变化[i]);
			cout << "位置量变化[i] x:" << 位置[i].x << " y:" << 位置[i].y << " z:" << 位置[i].z << endl;
		}
	}

}














float32 S_PBF解算器::f_解算密度(uint64 当前ID, const vec3* x, const float32* 质量, const float32 密度0, float32* 密度) {
	float32 错误密度 = 0;

	/*密度[当前ID] = 质量[当前ID] * m_粒子数据.m_GPU数据.m_Weight0;
	if (当前ID == 1891) {
		x[当前ID];
	}
	uvec3 coord = f_浮点坐标转索引坐标(x[当前ID], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移, m_邻接.m_Data.m_域属性.m_单元属性.m_单元大小);
	coord = vec_Max(coord, { 1,1,1 });
	for (int8 oz = -1; oz < 2; ++oz) {
		for (int8 oy = -1; oy < 2; ++oy) {
			for (int8 ox = -1; ox < 2; ++ox) {
				uvec3 搜索坐标;
				搜索坐标.x = coord.x + ox;
				搜索坐标.y = coord.y + oy;
				搜索坐标.z = coord.z + oz;


				S_DomUnitIndex 单元索引 = m_邻接.f_get单元粒子索引(搜索坐标);

				for (uint8 j = 0; j < 单元索引.num; ++j) {
					uint64 k = 单元索引.index[j];

					密度[当前ID] += f_权重(x[当前ID] - x[k]);
				}
			}
		}
	}
	//for (uint64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
	//	密度[当前ID] += f_权重(x[当前ID] - x[i]);
	//}

	错误密度 = DEF_Max(密度[当前ID], 密度0) - 密度0;*/
	return 错误密度;
}

float32 S_PBF解算器::f_计算拉格朗日乘数(uint64 当前ID, const vec3* x, const float32* 质量, float32 密度, float32 密度0) {
	const static float32 eps = 1.0e-6;
	float32 lambda = 0.0f;

	const float32 约束 = DEF_Max(密度 / 密度0 - 1.0, 0.0);

	if (约束 != 0.0f) {
		float32 梯度和 = 0.0f;
		vec3 邻接梯度 = {};


		/*uvec3 coord = f_浮点坐标转索引坐标(x[当前ID], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移, m_邻接.m_Data.m_域属性.m_单元属性.m_单元大小);
		coord = vec_Max(coord, { 1,1,1 });
		for (int8 oz = -1; oz < 2; ++oz) {
			for (int8 oy = -1; oy < 2; ++oy) {
				for (int8 ox = -1; ox < 2; ++ox) {
					uvec3 搜索坐标;
					搜索坐标.x = coord.x + ox;
					搜索坐标.y = coord.y + oy;
					搜索坐标.z = coord.z + oz;


					S_DomUnitIndex 单元索引 = m_邻接.f_get单元粒子索引(搜索坐标);

					for (uint8 j = 0; j < 单元索引.num; ++j) {
						uint64 k = 单元索引.index[j];

						vec3 梯度 = f_梯度(x[当前ID] - x[k]) * (-质量[k] / 密度0);
						梯度和 += vec_平方范数(梯度);
						vec_sub(&邻接梯度, 梯度);
					}
				}
			}
		}*/

		/*for (uint64 i = 0; i < m_粒子数据.m_动态粒子总数; ++i) {
			vec3 梯度 = f_梯度(x[当前ID] - x[i]) * (质量[i] / 密度0);
			梯度和 += vec_平方范数(梯度);
			vec_sub(&邻接梯度, 梯度);
		}*/

		梯度和 += vec_平方范数(邻接梯度);
		lambda = -约束 / (梯度和 + eps);
	}
	
	return lambda;
}




vec3 S_PBF解算器::f_解算密度约束(uint64 当前ID, const vec3* x, const float32* 质量, const float32* 拉格朗日, float32 密度0) {
	vec3 位移步进 = {};

	/*uvec3 coord = f_浮点坐标转索引坐标(x[当前ID], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移, m_邻接.m_Data.m_域属性.m_单元属性.m_单元大小);
	coord.x = DEF_Max(coord.x, 1);
	coord.y = DEF_Max(coord.y, 1);
	coord.z = DEF_Max(coord.z, 1);
	if (当前ID == 1891) {
		cout << "C x:" << coord.x << " y:" << coord.y << " z:" << coord.z << endl;
	}
	for (int8 oz = -1; oz < 2; ++oz) {
		for (int8 oy = -1; oy < 2; ++oy) {
			for (int8 ox = -1; ox < 2; ++ox) {
				uvec3 搜索坐标;
				搜索坐标.x = coord.x + ox;
				搜索坐标.y = coord.y + oy;
				搜索坐标.z = coord.z + oz;


				S_DomUnitIndex 单元索引 = m_邻接.f_get单元粒子索引(搜索坐标);

				for (uint8 j = 0; j < 单元索引.num; ++j) {
					uint64 k = 单元索引.index[j];

					vec3 邻接梯度 = f_梯度(x[当前ID] - x[k]) * (-质量[当前ID] / 密度0);
					邻接梯度 = 邻接梯度 * (拉格朗日[当前ID] + 拉格朗日[k]);
					vec_sub(&位移步进, 邻接梯度);

					if (fabs(位移步进.x) >= 6 || fabs(位移步进.y) >= 6 || fabs(位移步进.z) >= 6) {
						cout << "C x:" << 位移步进.x << " y:" << 位移步进.y << " z:" << 位移步进.z << endl;
					}
				}
			}
		}
	}*/

	return 位移步进;
}

void S_PBF解算器::f_更新速度矢量() {

	vec3* 速度 = m_粒子数据.m_速度.data();
	auto* 重量 = m_粒子数据.m_质量.data();
	vec3* 位置 = m_粒子数据.m_位置.data();
	auto* old位置 = m_粒子数据.m_上一帧位置.data();


	if (m_使用GPU计算) {

	}
	else {

	}

#pragma omp parallel for
	for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
		if (重量[i] != 0.0f) {
			//速度[i] = (位置[i] - old位置[i]) * (1.0f / m_更新步长);
			速度[i] = (位置[i] - old位置[i]) * 0.015;
			//速度[i] = (位置[i] - old位置[i]) * (m_粒子数据.m_粒子半径 / m_更新步长);
			//速度[i] = (位置[i] - old位置[i]) * (m_更新步长 / m_粒子数据.m_粒子半径);
			//速度[i] = (位置[i] - old位置[i]) * (m_更新步长);
			//std::cout << "m_更新步长 = " << m_更新步长 << " x:" << 位置[i].x << " y:" << 位置[i].y << " :z" << 位置[i].z << std::endl;
		}
	}
}

void S_PBF解算器::f_解算粘度() {
	vec3* 位置 = m_粒子数据.m_位置.data();
	vec3* 速度 = m_粒子数据.m_速度.data();
	auto* 密度 = m_粒子数据.m_密度.data();
	auto* 质量 = m_粒子数据.m_质量.data();
	
	float32 粘度 = m_粒子数据.m_粘度;

	if (m_使用GPU计算) {

	}
	else {

	}


//#pragma omp parallel for
	/*for (int64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
		uvec3 coord = f_浮点坐标转索引坐标(位置[i], m_邻接.m_Data.m_域属性.m_单元属性.m_域世界偏移, m_邻接.m_Data.m_域属性.m_单元属性.m_单元大小);
		//m_邻接.f_get单元粒子索引(位置[i]);
		coord = vec_Max(coord, { 1,1,1 });
		for (int8 oz = -1; oz < 2; ++oz) {
			for (int8 oy = -1; oy < 2; ++oy) {
				for (int8 ox = -1; ox < 2; ++ox) {
					uvec3 搜索坐标;
					搜索坐标.x = coord.x + ox;
					搜索坐标.y = coord.y + oy;
					搜索坐标.z = coord.z + oz;

					//uint64 index = m_邻接.f_从坐标取索引(搜索坐标);
					S_DomUnitIndex 单元索引 = m_邻接.f_get单元粒子索引(搜索坐标);

					for (uint8 j = 0; j < 单元索引.num; ++j) {
						uint64 k = 单元索引.index[j];

						vec3 v = ((速度[i] - 速度[k]) * (粘度 * (质量[k] / 密度[k]))) * f_权重(位置[i] - 位置[k]);
						vec_sub(&速度[i], v);
					}

				}
			}
		}
	};*/

}


//void S_PBF解算器::f_构建静态网格碰撞划分() {}

void S_PBF解算器::f_计算静态网格碰撞() {
	if (!m_粒子数据.m_静态三角网格顶点.size()) return;


	auto* v1 = m_粒子数据.m_上一帧位置.data();
	auto* v2 = m_粒子数据.m_位置.data();
	uint32 num = m_粒子数据.m_位置.size();

	vec3* 顶点 = m_粒子数据.m_静态三角网格顶点.data();
	uvec3* 索引 = m_粒子数据.m_静态三角网格索引.data();

#pragma omp parallel for
	for (int64 i = 0; i < num; ++i) {
		vec3 dir = (v2[i] - v1[i]);
		//float32 线段长度 = vec_len(dir);
		float32 线段长度 = vec_normalize(&dir);

		std::vector<S_八叉树子空间*> pathTree;
		f_空间划分_线段查子空间(m_粒子数据.m_静态网格空间划分, { v1[i], v2[i] }, pathTree);

		//std::vector<vec3> 相交点; 相交点.reserve(64);
		vec3 最近点 = v1[i];
		float32	最近点距离 = 1000000000;
		//uint32	最近点射线索引 = 0;
		uvec3	最近点三角形索引 = {};
		bool	是否有交点 = false;


		uint32 PathNum = pathTree.size();
		auto* p_path = pathTree.data();
		for (uint32 单元索引 = 0; 单元索引 < PathNum; ++单元索引) {
			uint32 单元数据num = ((std::vector<uint32>*)(p_path[单元索引]->m_数据))->size();
			uint32* indexs = ((std::vector<uint32>*)(p_path[单元索引]->m_数据))->data();

			for (uint32 k = 0; k < 单元数据num; ++k) {
				uvec3 index = 索引[indexs[k]];

				vec3 交点 = {};
				vec3 edge1 = { 0.0f, 0.0f, 0.0f };
				vec3 edge2 = { 0.0f, 0.0f, 0.0f };
				edge1 = vec_sub(顶点[index.y], 顶点[index.x]);
				edge2 = vec_sub(顶点[index.z], 顶点[index.x]);

				int32 是否相交 = f_graph_射线三角形相交(顶点[index.x], edge1, edge2, v1[i], dir, &交点);
				if (是否相交 > 0) {
					float32 dis = vec_len(v1[i] - 交点);
					if (dis < 最近点距离 && dis <= 线段长度) {
						//最近点射线索引 = i;
						//最近点三角形索引 = indexs[k];
						最近点三角形索引 = index;
						最近点 = 交点;

						是否有交点 = true;
					}
				}
			}
		}

		//f_graph_距离最近坐标(v1[i], 相交点.data(), 相交点.size());
		if (是否有交点) {
			//uvec3 index = 索引[最近点三角形索引];
			vec3 edge1 = vec_sub(顶点[最近点三角形索引.y], 顶点[最近点三角形索引.x]);
			vec3 edge2 = vec_sub(顶点[最近点三角形索引.z], 顶点[最近点三角形索引.x]);
			vec3 N = vec_normalize2(vec_cross(vec_normalize2(edge1), vec_normalize2(edge2)));
			vec3 反弹向量 = vec_Reflect(dir, N);

			float32 dis = vec_len(v2[i] - 最近点);
			反弹向量 = 反弹向量 * (dis + 0.00001f);
			v2[i] = 最近点 + 反弹向量;
			//v2[i] = 最近点;
		}
		

	}
}




float32 S_PBF解算器::f_更新解算步长(float32 最小步长, float32 最大步长) {
	vec3* 加速度 = m_粒子数据.m_加速度.data();
	vec3* 速度 = m_粒子数据.m_速度.data();


	//m_粒子数据.m_粒子半径;
	float32 h = m_更新步长;
	float32 系数 = 1.0;

	float32 最大速度 = 0.1;
	float32 直径 = m_粒子数据.m_粒子半径 * 2;
	for (uint64 i = 0; i < m_粒子数据.m_流体粒子总数; ++i) {
		float32 平方范数 = vec_平方范数(速度[i] + 加速度[i] * h);
		if (平方范数 > 最大速度) {
			最大速度 = 平方范数;
		}
	}
	h = 系数 * 0.4 * (直径 * sqrt(最大速度));
	h = DEF_Min(h, 最大步长);
	h = DEF_Max(h, 最小步长);

	m_更新步长 = h;

	return m_更新步长;
}

float32 S_PBF解算器::f_权重(vec3 r) {
	float32 rl = vec_normalize(&r);
	const float32 q = rl / m_粒子数据.m_粒子光滑半径;


	float32 权重 = 0.0f;
	if (q <= 1.0 && q > 0.0000001) {
		if (q <= 0.5) {
			const float32 q2 = q * q;
			const float32 q3 = q2 * q;
			权重 = m_粒子数据.m_GPU数据.m_K * (6.0f * q3 - 6.0f * q2 + 1.0f);
		}
		else {
			权重 = m_粒子数据.m_GPU数据.m_K * (2.0f * pow(1.0 - q, 3));
		}
	}
	
	return 权重;
}

vec3 S_PBF解算器::f_梯度(vec3 r) {
	vec3 梯度 = {};

	float32 rl = vec_normalize(&r);
	const float32 q = rl / m_粒子数据.m_粒子光滑半径;
	if (q <= 1.0f) {
		if (rl > 0.0000006f) {
			vec3 q梯度 = r * (1.0f / (rl * m_粒子数据.m_粒子光滑半径));
			if (q <= 0.5) {
				梯度 = q梯度 * (m_粒子数据.m_GPU数据.m_L * q * (3.0f * q - 2.0f) );
			}
			else {
				const float32 系数 = 1.0 - q;
				梯度 = q梯度 * (m_粒子数据.m_GPU数据.m_L * (-系数 * 系数));
			}
		}
	}

	return 梯度;
}

S_cuAdjacencyIndex S_PBF解算器::f_getGPU数据() {
	S_cuAdjacencyIndex cuAdjacency = {};

	/*cuAdjacency.m_Indexs = m_邻接.m_Data.m_邻接索引_Value;
	cuAdjacency.m_VoxelIndexs = m_邻接.m_Data.m_体素索引;
	cuAdjacency.m_IndexNum = m_邻接.m_Data.m_单元内邻接数量;

	cuAdjacency.m_粒子半径 = m_粒子数据.m_粒子半径;
	cuAdjacency.m_粒子光滑半径 = m_粒子数据.m_粒子光滑半径;
	cuAdjacency.m_Poly6 = 315.0f / (64.0f * 3.141592f * pow(m_粒子数据.m_粒子光滑半径, 9));
	cuAdjacency.m_Spiky = -45.0f / (3.141592f * pow(m_粒子数据.m_粒子光滑半径, 6));
	//cuAdjacency.m_Viscosity = 45.0f / (3.141592f * pow(m_粒子数据.m_粒子光滑半径, 6));
	cuAdjacency.m_Viscosity = m_粒子数据.m_粘度;
	//cuAdjacency.m_Weight0
	cuAdjacency.m_K = m_粒子数据.m_GPU数据.m_K;
	cuAdjacency.m_L = m_粒子数据.m_GPU数据.m_L;*/
	return cuAdjacency;
}

S_域区* S_PBF解算器::f_get域属性() {
	//return &m_邻接.m_Data.m_域属性;
	return nullptr;
}










